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The feedback stabilization problem for ensembles of coupled spin 1/2 systems is discussed from 
a control theoretic perspective. The noninvasive nature of the bulk measurement allows for a fully 
unitary and deterministic closed loop. The Lyapunov-based feedback design presented does not 
require spins that are selectively addressable. With this method, it is possible to obtain control 
inputs also for difficult tasks, like suppressing undesired couplings in identical spin systems. 
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I. INTRODUCTION 

NMR spectroscopy deals with the manipulation of nuclear spins of quantum ensembles, see These systems 

exhibit most of the essential features of quantum mechanical systems, like the state space of tensorial type (providing 
exponential growth of the degrees of freedom available) and natural coupling mechanisms between spins, which 
guarantee the nonclassical nonlocality characteristic of quantum evolutions. For the purposes of state manipulation, 
over the last 40 years the field of NMR has developed an extremely versatile and universally accepted set of tools, 
in the form of sequences of electromagnetic pulses 0, 0, 0| . In terms of classical control theory, these would be 
classified as open loop control methods. From the perspective of control theory, open-loop methods are by far more 
complicated and less robust than closed-loop methods, in which one or more functions of the state variables are measured 
on-line and used to impose a behavoir to the system through a feedback algorithm, or simply to reject errors and 
disturbances deviating the state from its desired trajectory. As new directions of research like quantum information 
processing are pushing for more precise and efficient engineering of quantum states, the use of a quantum/classical 
interface for the purpose of feedback control of a quantum state is becoming a reasonable task in different settings 
[Til H?^ 21, 27. 28] . From this control perspective, NMR systems constitute a remarkable opportunity to device 
feedback control methods at a quantum level for a number of reasons: 

1. the model of the system is known in detail; 

2. its control mechanism is also very accurate; 

3. the measurement is classical, thus avoiding all complications due to the state collapse problem (weak or less) 
unavoidable in other quantum control contexts |24L ; 

4. the relaxation times are sufficiently long to make the real-time interface with a control device feasible. 

For all these reasons, a completely classical, unitary and deterministic feedback in the context of NMR systems is 
theoretically feasible, although major technical problems remain to be addressed, like the very low signal-to-noise ratio, 
the on-line extraction of measurements from coil magnetization in real time and in presence of rf excitation and the 
need to on-line reconstruct the density matrix from the collective measurements. 

The aim of a feedback synthesis is to render a desired state or state trajectory an attractor for the system regardless 
of the initial condition 32], and the main purpose of this work is to investigate in detail this feedback synthesis from 
a theoretical viewpoint. The main tools we shall use are Lyapunov-based feedback design for bilinear control systems, 
adapting a method known as Jurjevic-Quinn condition 0, [I3 to the case at hand. Standard references for basic 
material on feedback control are e.g. the books 0,|2^. It is worth mentioning that in the physics literature a related 
approach (using model-based feedback only for the purposes of attaining the open loop control corresponding to a 
certain task, see also Sec. IVI|) appears under the name "tracking control" .22i. .30.] . 

For the purpose of feedback synthesis, the system is formulated as a bilinear control system living on a compact 
homogeneous space. For the task of tracking a given orbit, a class of control Lyapunov functions is naturally defined 
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by the notion of distance induced by the real Euclidean structure with which the homogeneous space is endowed. This 
construction resembles closely the Jurdjevic-Quinn stabilization technique (see also 0, 0, |2§| for related 
material dealing with pure states only), although the computation of the largest invariance set via LaSalle principle 
is more complicated. This is a consequence of the nontrivial topology of the state space (for a single spin 1/2 it is a 
sphere §^), implying that no (smooth) feedback design can achieve global stabilization (see for a detailed analysis in 
this direction). At most one can achieve convergence out of a singular set of isolated, repulsive points. The emphasis 
on the exact knowledge of the singular locus is motivated by the fact that near a singularity the convergence can be 
very slow. Moreover, a local design is of limited practical interest in a quantum context. For multispin systems it is 
shown that the tensor product nature of the state space does not complicate exceedingly the feedback synthesis. On 
the contrary, the singular set of the control law can be computed explicitly thanks to this tensorial structure. 

Throughout the paper we consider only the case of spins that are not selectively excitable. This is clearly the most 
difficult case, as an rf field affects all spins and interacts with all couplings. A similar feedback synthesis for selective 
controls is much simpler (especially for what concerns the convergence analysis) and can be deduced by similar means. 

While the analysis is easier for the Ising Hamiltonian, all the results are valid for different types of interactions like 
Heisenberg or dipole-dipole. In particular, we show (Sec. IIV|I how it is possible to reject unwanted coupling terms, 
provided they are sufficiently slow compared to the residual nonlocal part of the Hamiltonian. 

As an example of a difficult task that can be solved in this way, we consider in Sec. I VII a system of 3 identical spins 
in which we want to cancel the interaction between first and third spins without altering the dipole-dipole coupling of 
the linear chain. This is a typical problem araising in solid state NMR and for which a pulse sequence is not known. 
The feedback scheme allows to compute the time course of the control field achieving an almost exact decoupling in 
the entire Hilbert space. 



II. FORCED LIOUVILLE EQUATION FOR SPIN 1/2 SYSTEMS: A BILINEAR CONTROL MODEL 

The Liouville-von Neumann equation for a density operator p is: 

p^-t[H,p]. (1) 

For a single spin 1/2 system, assume the Hamiltonian H is composed of a free part (the drift, often called the Zeeman 
Hamiltonian) and a forcing part (the control term). If Aj, j = 0, . . . ,3, are the (normalized: tr(AjAfe) = Sjk) Pauli 
matrices, the free part, which is due to a strong static magnetic field Bo, is conventionally aligned with the A3 axis 
and causes the spin ensemble to precess around A3. In the laboratory frame this is 

Hf,e = -"fBoXs 

where 7 is the gyromagnetic ratio. The quantity tOo = jBo is normally referred to as Larmor frequency. The control 
Hamiltonian originates from an electromagnetic field rotating in the (Ai, A2) plane at a frequency o^rf close or equal 
to the precession frequency lUo- In the laboratory frame, this corresponds to the Hamiltonian: 

-ffrf = -7-B1 (cos(wi.fi + 4>)Xi + sin(wi.ft + (j))X2) 

where is the phase of the field. Typically, oji — jBi is of the order of the tens to hunderds of kHz. The controllable 
parameters are the amplitude Bi, the frequency Wrf and the phase (f). The resulting motion in the laboratory frame 
is rather complicated to describe. It is convenient to express it in a rotating frame, i.e., a coordinate system rotating 
about the vertical axis A3 at the frequency Wrf. It can be obtained by means of a variation of constants formula. 
Up to a negligible term, the neat result is that the inertial Hamiltonian Hf i + i/if is replaced by the rotating frame 
Hamiltonian H{ + He given by 

Hi = -(cJo-Wrf)A3 =/l^A3 (2) 

He ~ — (cos0Ai + sin(/)A2) . (3) 

Fixing (f) means fixing the axis at which the control acts. Hereafter we assume (j) = 0. Hence becomes simply: 

He = uXi (4) 

with u = —uji as our real valued control u G C°°(R). When uj^i ~ ojo, the Hamiltonian is driftlcss and the resulting 
motion is called nutation around the axis determined by (j). 



3 



Consider a two spin 1/2 system in an external field. To represent the Hamiltonian, we make use of the prodout 
operators basis given by Kjj^ = Aj (8) A^, j, /c = 0, . . . , 3. Calling 7q, and 7^ the gyromagnetic ratios of the two spins, 
the r.f. Hamiltonian now becomes 

-ffrf = -Bi (C0s(cji.ft + 0) (7aAio + 7/3A01) 

+ sin(cjrfi + (j)) (7qA2o + 7/3^02)) 

or, in the rotating frame described above and with = 

He = -Bi (7aAio + 7/3A01) . 

In the same rotating frame, the free part of the Hamiltonian H{ is 

Ht = h"^Ao3 + h^^Aso + /i''A33 (5) 

where h'^^ and h^^ are the differences between the Larmor frequencies of each spin, call them LUo,a and Wq^^, and the 
carrier frequency Wrf, h^^ = —{LJo,a — i^ii), h"^ — — (wo./s — i^rf), and h^^ represents the so-called J (or scalar) coupling, 
with typical values of hundreds of Hz. 

If the spins are homonuclear, 7^ = 7/3, then ujo.a and ujo.p (and thus h'^^ and h'^'^) differ only because of the chemical 
shift, typically of the order of a few kHz (compared to the MHz order of both Wo,a and Wo,/?)- If instead we have 
heteronuclear species, then this difference may be of the order of several MHz. 

We shall distinguish between nonselective and selective Hamiltonians. The main difference between the two cases 
is in the structure of the control Hamiltonian He. In the first case, it is assumed that the different spins all experience 
the effect of the same control field 

iJc = uHc,ns = (Aoi + Aio) (6) 

where u = — 7qBi = —jpBi. This model is suited for homonuclear species with LOo,a and iOo^fs not sufficiently 
separated. When instead tOo^a and cJq.^ are well separated, like in homonuclear species with consistent chemical 
shift or in heteronuclear species, it is a good approximation to consider the spins as selectively excitable |33|. This 
is achieved by considering 2 different rf fields tuned around the two frequencies ujo,a and uJo,i3 call them Wrf^a and 
<^rf,/3- The separation of Wo,a and Wo,/? implies that the cross-talk with the off-resonant spin is negligible. By applying 
the double coordinate change, one gets still the free Hamiltonian |5j, but now with = —{uJo,a — ^ii.a) and 
/i"^ = —(ujo,0 — (^ri,f3)- Choosing Wrf.a = ^o,a and uji-[^p — u}o,(3, the two local precessions disappear from Q. The 
coupling term is unchanged, as A33 commutes with A03 and A30. The control Hamiltonian, instead, becomes 

He = uoiAoi -1- wioAio (7) 

where mqi = -IpBi^p and uio = ~laBi^a- 

Models like (jSJ with only a vertical coupling are often referred to as Ising Hamiltonians. In particular, notice that 
since Hi is diagonal, 

rft 03 ^ _,. /,33 1 
Hf - _ j^ZO _ ;j33 

when h^'^ — h?^ jsj] the unforced system has at least one degenerate eigenvalue of multiplicity 2, regardless of the 
value of h^^ . This may result in a loss of controllability and complicates also the convergence in the closed-loop 
system. The diagonal elements of Hi have the meaning of energy levels of the (unperturbed) system. From 4] , since 
Graph(iJc,ns) is connected, as soon as Hf is -ffc,ns-strongly regular, i.e., has energy levels that are nondegenerate and 
transition frequencies all different in correspondence to the nonzero elements of -ffc.ns, then the system is controllable, 
see Theorem 3 of 0. 

Lemma 1 Consider the system Hi is Hc^ns- strongly regular if hP'^ 7^ li?^ , li?^ 7^ and Ir''^ 7^ ±(/i"'^ — h^^)/2. 

Proof. The graph of the control Hamiltonian iJc,ns habilitates the following 4 (nonoriented) transitions: 1^2, 
1 ^ 3, 2 ^ 4, 3 ^ 4. Computing the energy differences in terms of the h'~'^, h^'~' and h^^, lack of degenerate 
transitions corresponds to the three inequalities stated above. □ 

The condition of the Lemma is a generic condition, satisfied almost always. Furthermore, it is sufficient but not 
necessary for controllability. 
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Corollary 1 The system is controllable if ^ , h^^ ^ and h^^ ^ ±{h°^ - h^°)/2. 

For n spin 1/2, in a rotating frame of frequency oj^i, the Ising Hamiltonian of a linear spin chain is still composed of 
a drift part containing the Larmor precessions (relative to u^i as in plus the J couplings between adjacent spins 

and of a forcing term along the Ai axis of each spin which, in the nonselective case, is 

He = uHc,ns = U (A0...OI + . . . + A10...0) ■ 

Rather than working with the complex matrix ODE (whose integral flow corresponds to a conjugation action) 
we prefer to use the real "one-sided" linear action on vectors (or tensors) of expectation values along a complete 
orthonormal set of observables. For a single spin 1/2, it is well-known that one can write the density operator p as 
a real vector of expectation values along the Pauli matrices Xj, j = 0, . . . ,3: p = p' Xj = g ■ X, using the summation 
convention. In terms of g, and for the Hamiltonian in Q and Q one obtains the Bloch equations: 

g ^ -i (h^adx^ + uadx^) g. (8) 

The notation "ad" in ^ originates from the notion of adjoint representation, and the adjoint operators ad^^ stand 

for matrices of structure constants with respect to the su(2) basis given by the —iXj: adA^ A^ = [A^ , Afe] = Ym=o '^\k^i- 
In general, the adjoint representation of a semisimple Lie algebra is a real isomorphic matrix representation of the 
algebra. For su(2) for example, we have ad5u(2) = so (3). The most important feature of the adjoint representation 
in our case is that it provides a linear representation of one-parameter groups of automorphisms of the algebra. This 
enables us to formulate the control problem in terms of standard real bilinear control systems also for multispin 
systems. In fact, for 2 or more spin 1/2 densities, a parametrization similar to the Bloch vector yields a tensor, called 
the Stokes tensor and also denoted by g: if p G H®^, p = g^^Ajk — g - A where p''^ = tr (pAjk) are expectation values 
capturing all 15 degrees of freedom of p along the complete set of observables A = {Ajfe, j,k — 0, . . . ,3}, g £ S C S^^, 
where the real set S is very difficult to describe explicitly. All details are given in [J- Calling 

adAjfc = ^ (adAj «) aadA^ -I- aadA^. adA^) 



the real skew-symmetric operators obtained from the adA^ above and the "antiadjoint" operators aadA^ (which have 

Afe} = ELo4fc' 



a similar meaning, only involving the "symmetric" structure constants aadA^Afe = {Xj, Afe} = X]f=o ^U-^Oi then we 



obtain the following adjoint representation of the Liouville equation 

^ = (ad/fj. + wadif^ „J £> (9) 

or, in components, 

pP"^ =-i (/i°3adAo3 + /i3"adA3„ + /i^^adA33)j 
-lu (adAoi +adAiJf^p'™. 



(10) 



By writing p^'^ as a 16-vector and expanding the tensor products, a bilinear control system with drift and control 
vector fields that are 16 x 16 matrices is obtained. The expression is similar for the selective case ©-0- 

Obviously, from g2s = Lie{— iAjfe, j,k = 0, . . . ,3} — su(2) 0su(2) Usu(2) (8)su(2), one gets for the adjoint represen- 
tation adg2^ — Lie{— iadAjfc , j,k = 0, . . . ,3} = so(3) © so(3) U so(3) (g) so(3). 

In the case u = the Ising model is completely integrable because 

[adAo3 , adA3o] = [adAo3 , adA33] = [adAgo , adA33] = 0- (11) 

Hence 

gP^it) = (exp {-ith°^adA,,) exp {- it h^'^ ad a,,) exp {-ith^^adA,,)Y^^^ ^''"(0) (12) 

which shows that the evolution is a "superposition" of 3 different periodicities, two local Tp^ ~ and r^^ — 
and one nonlocal Tp — The purity of the two reduced densities Pait) — tr^(p) and Pf3{t) = tra{p) is varying with 
time because of the coupling. 

The generalization to n spin 1/2 is completely analogous: the 2" x 2" density matrix p can be described by an 
n-index tensor p = gf^^ "^" Aj^ .. = ^ • A, each index ranging in 0, . . . , 3, Aj^ .. = Aj^ (g) . . . (g) Xj^ . The corresponding 
ODE is still given by an equation like 0. 
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III. STATE FEEDBACK STABILIZATION OF SPIN-1/2 SYSTEMS 

In this Section, we are only interested in full state feedback. The topology of the manifolds discussed in this work 
(spheres and compact homogeneous spaces obtained by taking the "envelope" of tensor products of "affine" spheres) 
forbids to have globally converging smooth algorithms. For example, for the Bloch sphere §^ there does not exist 
smooth positive definite functions with less than two points having zero derivative. In control theory, the functions 
having such minimal number of zeros are sometimes referred to as Morse functions |2Cll |. Hence in the simplest case 
of a single spin 1/2, the Lyapunov-based design will always be characterized by the presence of at least a spurious 
equilibrium point, which can, however, be rendered repulsive. 



A. State feedback stabilization for a single spin-1/2 system 

Because of the drift term, for q different from the north/south poles of S'^ the unforced system q — —ih^adx^Q 
does not have an equilibrium point but keeps moving on an orbit in the {gi, g2) plane corresponding to g^{t) = g^{0). 
We shall not try to cancel such precession motion, but instead try to stabilize the state to a desired orbit. Consider 
the system ^ from a given initial condition ^(0). If is the reference state, describe the orbit to track (with the 
obvious prerequisite ||£»d|| = ||^|| = r) by means of a ODE like ||SJ| but without forcing terms. Calling Qd{t) the desired 
reference state, 

Qd ^ -ihl&dx^Qd (13) 
means that g\ and g^ evolve on a circle while g^{t) — g'^{0) is the fixed value that characterizes the orbit. 

Proposition 1 The system ((SJ with the feedback law 

u = k{{Qd, -iadxiQ)), (14) 

where k G M'*', is tracking the reference orbit Qd{t) given by lO with h\ = h^, in an asymptotically stable manner, 
for all q{Q) £ with the exception of q{0) — — ^d(O) and of q{0), Qd{Q) such that g'^{0) — g'^ — 0. 

Proof. 

For density operators, a natural choice of distance is given by the so-called Hilbert-Schmidt (or trace) norm: tr (p^) . 
In terms of Bloch vectors ^d, Q, this induces the following distance between ^d and g, see e.g. .SlJ: 

d{gd, g) = ll^dlP - ({gd, g)) = llsdlP - gjg- (15) 

Consider as candidate Lyapunov function the distance (|15|l : V — d(gd,g)- Clearly V ^ and V — only when 
g,i = g. Since = 

v = -{{gd,g))-{{gd,e)) 

= -{{-ih^a,dx.,gd, g)) - {{gd, ~i (^^adA,, + uadAi) g)) (16) 
= iugjadx^g, 
because iadx^ is skew-symmetric. Inserting (I14II : 

V = -k{ig^a.dx,gf ^0. 

The system (0 with the feedback (I14II is time- varying and it is not possible to get rid completely of the time dependence 
by considering the error system e = ^d ^ 

e = — ladfffC — iuadn^ „s(e — gd)- 

However, from 1)16(1 . V is time-indepedent thanks to the cancellations and so must be V (its time-dependence is only 
apparent). Hence we can use LaSalle invariance principle for autonomous systems and the positive limit set of g{t) 
is still the largest invariant set, call it £, confined to Af = {g such that V = 0} (and corresponding to u = 0). To 
compute £, following the same idea of the proof of Theorem 2 of |0, in Af it must also be ^ = 0. Explicitely: 

^ = -ik{g^a,dx,g+ gladx.g) (17) 
= kh^g^[—iad\-^, —iad\^]g = V2kh^igjadx2g- 
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Notice, however, that u = ^ = yield bihnear forms as opposed to the quadratic forms of the original proof of 
[Tsi] : that in addition we have the constraint of ||^(i)|| — const ^ to deal with; and that the Lie algebra involved 
is composed of only skew-symmetric matrices, which applied to a point yields (out of the singularities) the tangent 
plane to the sphere §^ (not K^). This makes the condition of nonglobal. For example both bilinear forms l|14|) 
and H17I) are identically zero on the great circles g'^ — 
then in A/" we have 



0, regardless of the values of g'^, g^ , j = 1,2. If g'^ ^ 



6d 



J = 1,2. 



(18) 



and the closed loop system confined to J\f is given by the drift alone: 



= V2h^ 








• 

-9l 




(19) 



The presence of spurious equilibria in M is equivalent to the feasibility of 1)19(1 with g^ ^ g\. In A/", expanding the 



isospectral constraint ||^d|| = 



one gets: 

+ (A? 



^ {(sir + + [or = 10 {{sir + i^ir + i^i?) 



(20) 



The first and last expression of H20I) are equal if and only if j^jjr = 1, i.e., g^ = ^Qd- While g^ — +g'^ is already on 

the equilibrium since 118|l implies that ^ = ^d, the case g^ = —g'^ lead to a singular point of the control law given 
by Q{t) = — ^d(^) (again from t^l8\i ). Such a singularity corresponds to the antipodal point to the current desired 
position and is an isolated unstable equilibrium point. Hence, whenever 7^ 0, £ = {Q{t) — iffd(0} ^-^^ ^^'^ closed 
loop system almost globally tracks the desired orbit in an asymptotically stable manner. Since V is almost always 
decreasing along the trajectories of the closed loop system and the condition g^ does not depend on time, it is 
possible to recompute £ in terms of the initial conditions: 



u{t) = - i/c£.d (0)e"''''''*^^adAie-"'*'^^^3£,(0) 

^0 



sm{V2th^) cos{V2th^) 

=V2k (cos{V2th') {gl{0)g\0) - gl{0)g'{0)) 



=V2fceI(0) 





sm{y/2th^) 
coslV^th^) 




e(o) 



(21) 



,iV2th^) {gliO)g\0) ~ g\{0)g'iO))) = 



leads to conditions equivalent to p8|) and H20|l in terms of q{0) and £>d(0). Hence Q{t) = —Qd{t) if and only if 
^(0) = — 5d(0), i.e., the system is initialized in the antipodal point with respect to the initial condition ^d(O) of the 
desired trajectory. 

□ 



Remark 1 The closed loop system is a nonlinear Bloch equation very similar to those obtained by considering 
radiation damping effects 

Remark 2 A periodic solution of an autonomous system can never be a global attractor of a compact set (like the 
sector of the sphere delimited by ^ and not containing — ^d(O)). The time- varying formulation allows to bypass 
this topological obstruction. 

Remark 3 The exact cancellation of the drift in (|16|) is crucial for the proof of stability. If ^ h^, in fact, (|16|l is 
not homogeneous in u and the set of singular points in M is larger. 
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Remark 4 If rather than tracking a given orbit (a full state stabilization problem) we are interested only in the 
orbital asymptotic stabilization of the invariant set — const, then the problem becomes unidimcnsional. In this 
case, simple distance-like feedback laws like 

u^kg\gl-g'') (22) 

readily provide a solution. 

In Figg. [T151 simulations of the closed loop system with the controller (|14|l are shown. In Fig. O the instability of 
the antipodal point is shown: while q{Q) = —Qd{0) implies the state (dashed line) is not converging to Qd{t) (dotted 
line), a small perturbation in ^(0) is enough to make Q{t) (solid line) converging to Qd{t). 




FIG. 1: Closed-loop trajectory on the Bloch sphere for the controller 11411 . 



FIG. 2: The components of the Bloch vector of Qdit) (dotted line) and Q{t) (solid line) for the same data as in Fig.0 



FIG. 3: The closed-loop trajectories of the system with intial state (dashed line) "antipodal" to the desired state (dotted line) 
Qd{0) = — e(O) and from the same antipodal initial state plus a small perturbation (solid line). 



When the rf field rotates at the Larmor frequency, the Hamiltonian of the system is driftless as = ujo — ujif = 0. 
In this case, in the rotating frame we have a nutation motion, i.e., just a rotation around the Ai axis: 

Q = —iuadx^Q. (23) 
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Assuming Qd{t) = Qd{0), then both the feedback laws (|22(l and H14|l still hold and give a time-independent closed loop 
system. As a consequence, almost global stabilization is no longer achievable on a compact and in fact both schemes 
correspond to partial state stabilization schemes, since in H23|) ^ remains critically stable in closed loop. 

B. State feedback for 2 spin-i systems 

Assuming the entire two spin 1/2 state tensor g is available on-line, we want now to generalize the feedback law of 
Section lill Al to the ensembles of weakly coupled spin 1/2 described above. Only the nonselective case will be treated. 
Clearly, the selective Hamiltonian of eq. Q) has more control authority available, hence its stabilization is an easier 
task than the nonselective case. As before, let Qa, Qp and Pa^, Qp^ be the reduced densities respectively of g and 
^d- An obvious prerequisite for stabilizability is that g and both belong to S, which implies for example that 
II ^11 = ll^dll- We shall also assume that the intial condition is a product state ^(0) = gaiO) ® Qpi^)- 

The proof of Propositionnwas essentially relying on the Jurdjevic-Quinn sufficient condition for stabilizability [T^ . 
opportunely modified in order to deal with skew-symmetric infinitesimal generators. The key tool is the so-called 
"ad-condition", i.e., a particular type of Lie brackets often used for testing controllability, see Q. The following 
proposition shows that this condition is never satisfied by the system 

Lemma 2 For the system ^ 

spanj— iad/ff, —iadn^^^, [—iadHf, — ?ad/f^ ^J, . . . , [— iad/ff , ■ • ■ i ["^a-d/fj, — iad//^ ^^J . . .] } T^adg^^. (24) 

Proof. Given any A,B^ 02s, for the adjoint representation 

[ad^, ads] = ad[A^B], (25) 

with ad^jads G adg^^. Hence, (|24|) holds if and only if 

span{-ii/f, -ii/c,ns, [-iH{, -«i?c,ns], ■ • ■ , [-iHf, . . . , [~iH{, -ii?c,ns] ■■■]}¥' 32s- (26) 

Computing explicitly the first Lie brackets: 

[-iHf, -zHc.nJ = -i {h°^Ao2 - /i'^A32 + /i'°A2o - h^^A23) (27a) 
[-iHf, [-iHu -zi/cns]] = -I {{{h'' f - (h'^'f) Aoi + {{h'^^f - (h^y) Aio) (27b) 

Comparing JB)) and (|27b|) it is clear what the pattern of iterated Lie brackets will be, since the basis elements involved 
are the same. In particular, the basis directions — iAn, — zAi2, — iA2i and — iA22 are never touched by such sequences 
of commutators. Hence the claim (|26f) or, equivalcntly, H24I) . □ 

As for the single spin 1/2 case, we want to achieve asymptotically stable tracking of the following periodic orbit: 

£»d = - iad^, , , 

= - i {hfadAg^ + ft.d°adA3o + hfadA^^) gd, gd{0) G S 

by means of a single control input. 

Proposition 2 Whenever Hf is Hc^-as- strongly regular, the feedback 

u ^ k{{gd, -ia,dH,,^,g)) (29) 

with k G M^, asymptotically stabilizes the system to the time-varying reference state gdit) given by H28|) with 
Hf^ = Hf, for all g{Q) G S with the exception of g{0) such that (ect(O), gpi^)) = — (^^^(O), e/3d(0)) and of all pairs 
(e(0), ed(0)) having {gl, J ^ (q, 0) and {g^, g^J = (0, 0). 

Proof. As in the proof of Proposition ^ take as Lyapunov function the analogous of the distance (|15|l . V{t) = 
ll^dlP — {{Qd{t), Q{t)))- Again, when differentiating the drift disappears, 

V = ~((^d, Q))-{{Q<i, q)) (30a) 
= i (((adHt£>d, q)) + {{Qd, {adHf + uadn^^^Jg))) (30b) 
= -u{{gd, ~iadH, „,g)), (30c) 
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and V is made negative semidefinite by the choice of feedback (|29(l . Comphcations arise when checking the emptiness 
of the invariant set £ in Af = {g s.t. V = 0}. In fact, from Lemma[21 the Jurdjevic-Quinn condition never apphes 
to tensor product systems. Furthermore, since the reduced densities Qa{t) and Qi3{t) in Af have time- varying norm, 
neither the method used in the proof of Proposition ^ is directly apphcable. However, it can be apphed "strobo- 
scopically" i.e., at the time instants t = CTp, c € N, in which the nonlocal part of the dynamics disappears. From 



du 



where we have used (I25II . The explicit expression for ad[//^ „b,H!] follows applying (|25|l to H27a|> . From ()llfl and H12|l . 

in J\f one can write 



u{t) =-ik (eJ(0)e^*''''"'^^33e'*'''"^dA„3^^^^^g-»t/««^adA 



03 g-it/i^'^adAaa ^^0) 



(31) 



and similarly for ^ = 0. At t — cTp, g±«*'i^^adA33 _ j^^ hence (|31|) contains only local terms and the integration in g, 
gd resembles (PT|l : 



iV2k (eI(0)/4 ® e**''°'"^^3adAie-**''"'"^^3 ^(0) + gj{0)e'"' 
I 



■'ad; 



3 ad 



-ith3''adA, 



.74^(0) 



--2k 



el(0)ea(0)®el(0) 



V 



+eI(o) 









.0 sin(t/i03) cos(t/i03) 



-sin(i/i30) 

-cos(i/i30) 

.0 sin(i/i30) cos(i/i30) q 





sin(i/i03) 

COs(t/l03) 





(32) 



ea(0)®eI(0)e/3(0) =0. 



Eq. H32|) is enough to show that the structure of f at i = cr^ reflects essentially that of Proposition ^ In fact, 
the assumption of i7c,ns-strong regularity of imphes that (excluding the case (0)^^(0) — g^pJ^)gp(W) — 
corresponding to the antipodal case, but only for pure reduced densities), each of the two terms of (|32|) must be zero. 
Therefore in M we have the following invariant states (for t = CTp): 



1. ea(0) = — ffQd(O) ^^^d to £>/3(0) = — e/5d(0) ("local" antipodal case); 



2. 



.3—^3 _ 



= and/or q^^ = q^^^ — 0, regardless of the values of the other terms g"-', (P^ , g-j^ , j = 1,2 



(reduced densities of true and desired trajectories are horizontal great circles, for any degree of mixing, including 
the completely mixed case). 

All cases above correspond to u = and give the invariant set £ at t = cTp. Hence, for any t, £ will contain at most 
a subset of them. All are excluded by the assumptions of the Proposition. □ 

For a typical simulation, the entire 16-state reference tensor (dotted) and the tensor of the closed-loop system 
(solid) are shown in Fig. ^ and the two reduced densities in Fig. |S1 

Remark 5 Since — iad//^ „^ is local, the feedback action (|29|l preserves ||£»q|| and ||^^||. 

Remark 6 From a practical point of view, although out of the singular set convergence is guaranteed, the rate of 
convergence may be very slow for inital conditions near the singularities. Having a second nonselective control field 
at the same frequency but with a phase difference of 90° i.e., using 



iJc = wi(Aoi + Ai 



W2(Ao 



A, 



in place of ©, with the feedback law U2 = {{gd, — iad(Ao2+A2 
rate. 



o)g)) helps considerably in improving the convergence 
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FIG. 4: The 16-components of the tensor (dotted hnes) and g (sohd hne) for the orbit tracking problem. 




FIG. 5: Closed-loop reduced densities Qa{t) (left) and Q0{t) (right) corresponding to Fig. |1| 



In the selective case described by the Hamiltonian ((SJ and O, the two feedback laws for uqi and uio are of the same 
type as (|29ll . Tuning the two rf fields exactly on the two Larnior frequencies cuo.a and Wo./j, one gets h^^ = — 0. 
Therefore the proof of Proposition [3 does not apply straightforwardly. However, the idea behind the proof still 
holds true, namely that the two reduced closed-loop evolutions are nonidentical. The presence of two controls in the 
selective Hamiltonian Q emphasizes this distinguishability and simplifies the stabilizability argument considerably 
with respect to Proposition |21 

Another consequence of this "distinguishability" notion, as a generalization of the argument in the proof of Propo- 
sition[21 is that the same tracking scheme can be used for more complicated free Hamiltonians than In particular, 
any 

Hi = /i°3Ao3 + /i^^Aao + h^'Ajk, 3,k^0 (33) 
can be used in Proposition [3 provided that Hf^ = Hf. We state it as a Corollary. 

Corollary 2 If Hf given by H33|) is H^^ns- strongly regular, the feedback law (|29|l asymptotically stabilizes the system 
Q = —isLd(Hi+uHc „a)ff to the time-varying reference state Qdit) given by Q(\ — —indHiQd for all q{0) G S, except for 
the same singular set described in Proposition\^ 



C. n spin 1/2 case 



For the Ising Hamiltonian, if we start from a product state, label the spins as a, . . ., and call pa , . . Pv the 
corresponding reduced densities, then we can arrive at the same conclusion as Proposition |21 provided all energy 
levels are neither equal nor equispaced. The proof is omitted as it makes use of the same techniques used above, only 
the notation is more cumbersome. 
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Proposition 3 If Hf is Hc^na- strongly regular, the feedback 

u = k{{Qd, -zad_H-,,„3e)) (34) 

with k G K"*", asymptotically stabilizes the system q = — i(ad/jj +ua,dH^-^^)Q to the time-varying reference orbit 
Qd{t) given by = —i&dHiQd , for all q{Q) such that ||^(0)|| = ||gd(0)||, with the exception of the antipodal point 
(ea(0), . . . , ^1.(0)) = -(ead(O), . . . , QuA^)) and of all pairs (^(0), ed(0)) having {qI, glj = (0, 0), . . {gl, qIJ = 
(0, 0). 

IV. SUPPRESSING UNWANTED WEAK COUPLINGS 

The J-coupling used in the previous Sections is an indirect coupling mechanics, physically due to the electrons 
shared in the chemical bonds between the atoms. Apart from this coupling, there are other interaction mechanisms, 
due to the direct or electron-mediated interactions between the spins. 

Consider the nonselective two spin 1/2 Hamiltonian (|33() and 

Hi, = h°^Ao3 + ft'°A3o + h'd^jk, J, k^O. (35) 

Call Hs — Hi, — Hi the difference between the desired and the true Hamiltonian. We are interested in treating the 
extra terms Hg as disturbances and suppressing them by means of feedback. 

Assume the frequencies h']^ are of the same order of magnitude, call it 27r/Td. When the frequency of the disturbance 
Hs^ call it 271 /ts (~ /i]-*^), is much smaller than 27r/Td, then can be suppressed by the control action. 

Proposition 4 Assume Hi is Hc^ns- strongly regular and ts >> rj. Then there exists a tOrt and a sufficiently high 
gain k such that the system g — — iad(/fj_|_„^f^ ^^)Q with the feedback ()29|l can track the reference trajectory Qd{t) given 
by Qd — —iad/fj^ Qd and reject the disturbance Hg up to a bounded error. 

Proof. Since Hi, ^ Hi, in the proof of Proposition[21the derivative of the Lyapunov function is no longer homogeneous 
in the control. Inserting the feedback (^5)) . 

V = {{Qd, ~ia.dHsQ)) - k{{Qd, -iad//^ „^ (36) 

the first term is in general sign indefinite, implying negative semidefiniteness of V is not guaranteed. Consider 
the set Af — {g s. t. V — 0}. Since k{{Qd, —iadH^^^Q))'^ ^ 0, for an invariant set to belong to Af it must be 
((^d, —iSidHsQ)) ^ 0. As — iad//^ is skew-symmetric, the only solution is ((gd, —iS'^^H^Qf) = {{Qd, — «adffc„a^))^ = 0- 
This happens for example for the perfect tracking Q — Qd- When this is not verified, the steady state is replaced by 
a limit cycle which is in general stable but not asymptotically stable. 

The first term in (|36|l is slow, of small amplitude, of period ts (when u — 0) and of zero average. Provided k is 
sufficiently high, the second term has a fast dynamics with respect to the first one and a large amplitude. Hence, in 
the time scale Td, ((^d, —i&<^HsQ)) can be thought of as frozen, with the controlled term still acting as a damping. 
The cancellation of Hg cannot be perfect because there are regions of the state space where the controllable term 
of (|36|l vanishes, namely where any of the 4 reduced densities Qa,, Qp,, Qa, and approaches the Ai axis. Recall 
(Remark Isj that the amplitude of the feedback action only depends from the reduced densities. In absence of a local 
precession motion, the system leaves this non-convergence region only due to the coupling terms and it is not possible 
to guarantee a recovery from the disturbance-induced instability for all Hi,. However, since h^^ = — (wo.a — ^ri) and 
= —{lOo,i3 — Wrf), it is always possible to choose iVri so that 1/ts is small compared to h^^ and h^^. The effect of 
the local precessions is to steer the corresponding reduced dynamics, both the desired and the real ones, out of the 
uncontrollable alignment with the Ai axis. Therefore, since both local closed-loop dynamics evolve fast compared to 
Hg and so does Hi,, the displacement due to the drift term in (|36|l can be rejected in the fast time scale. Since the 
disturbance Hg is persistently exciting the system, no steady state is ever reached, but the error remains bounded. □ 

Two common coupling models often used in the literature are the Heisenberg interaction and the dipole-dipole 
interaction 0, 0| . For example, the dipole-dipole Hamiltonian is given by 

Hdd = -t^dd(Aii + A22 - 2A33) (37) 

and models the direct coupling between the magnetic moments in solid state NMR. As an application of Proposition^] 
we want to modify the response of Hdd so as to reproduce that of a J coupling. In the rotating frame Wrf, assume 
that the free Hamiltonian is 

Hi = /i"3Ao3 + /i'°A3o + /i'^A33 + /i"An + /I'^Ajz (38) 



12 



where Hs = /i^^An + h?'^K22, h}^ = K^"^ = = -Wdd and lOwdd ^ hP^ ^ h?^ ■ The control Hamiltonian is 

stih nonselective and given by ((HJl. A typical closed loop behavior for this choice is shown in Fig. Even after the 
offset due to the initial condition is recovered, the system does not quicky reach an unperturbed steady state due to 
the periodic excitation given be ■ The Lyapunov function V and its derivative V are shown in Fig. [71 Notice 




Vvws ^f^frntx "^-fiMfWi /VV/a- 



10 20 10 20 10 20 10 20 



FIG. 6; The 16-components of the tensor (dotted lines) and g (solid line) for the weak coupling suppression problem. 




5 10 15 20 25 



FIG. 7: Top: the Lyapunov function V. Bottom: V (solid line) and its 2 components (dotted and dashed lines), see l|36^ . 



further that Corollary |21 implies that the coupling included in Hf^ needs not be restricted to the vertical A33 direction. 
Only the weakness of the unwanted coupling H{ — with respect to and with respect to the local precession 
frequencies matters for the disturbance rejection. 



V. FEEDBACK FROM MEASURABLE QUANTITIES 

A typical NMR measurement apparatus can provide a collective magnetization measurement in the (Ai, A2) plane. 
We assume this measure process can cohexist with the excitation of the coil (actuator signal). While for a single 
spin ensemble this, together with \\q\\ = const, allows to easily recover the entire state vector, for two spin systems 
what is measurable depends on the nuclear species we are considering. In the homonuclear case, the measurement 
corresponds to the output vector y — [yi 2/2]"'", where yj{t) = g^^ {t) + g^'^{t). Since also the controls are nonselective, 
a feedback from y cannot stabilize the system. However, if we can, by means of a state observer 35], access the 
entire 3- vectors and g^'^ rather than just the sums of their (Ai, A2) components, then several control designs are 
possible. We will not investigate the state estimation issue in detail here, see e.g. [l3- We only notice that the 
system is easily seen to be observable as soon as it is controllable but that due to the complicated structure of 
S, linear state observers are inadequate (they do not preserve ||^||). In the heteronuclear case, the estimation of 
and Qp is an easier task than in the homonuclear one. In fact, for different spin species, different coils (tuned at the 
different, well-separated, precession frequencies Wo.a, ^o.p) can be used so that both pairs {g''^'^, p"^) and (g^", g^") can 
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be available simultaneously from direct measurements. Since in NMR experiments the initial condition q(0) is always 
known, the simplest way to recover the g'^^ and g^^ components in this case is obviously to numerically integrate the 
system. 

Assuming Qa{t) and Qf3{t) are available, then one could use Qa{t) (81 Qpit) in place of Q{t) in the feedback controller 
of Proposition 12 

u = k{{Qd, -~ia.dH,^„,Qa <^ Qis))- (39) 

The approximation of the true state Q{t) with the product state Qa{t) ® Q/^it) corresponds to disregarding at each 
time the correlation that is being built by the J-coupling h^^. From Remark |S1 this difference is to a large extent 
negligible. The results are plotted in Figg. I8I9I comparing with Figg. 141 5 1 notice how the transient is only slightly 
worsened. 




-1 1 1 I 1 I 1 I 

10 20 10 20 10 20 10 20 



FIG. 8: The 16-components of the tensor Qd (dotted lines) and g (solid line) for the orbit tracking problem with feedback from 
reduced density measurements. 




FIG. 9: Closed-loop reduced densities Qa{t) (left) and Q/iit) (right) corresponding to Fig. |Sl 



VI. DESIGN OF OPEN LOOP CONTROL VIA MODEL-BASED FEEDBACK 

The feedback schemes of the previous Sections have a straightforward application as open-loop model-based state 
steering methods, see Fig. ^] In this setting, the system and its feedback controller are simulated on a computer. 
The time- varying control signal u{t) obtained in this way can then be used in an open-loop fashion to control the 
true system in the lab. This strategy is particularly significant for very difficult control tasks, like canceling unwanted 
couplings in spin systems not admitting selective controls, like is often the case in solid state NMR. Consider the 
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feedback law 
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FIG. 10; Block diagram for the "model-based feedback" control design. 




FIG. 11: Three spin configuration. 



3 identical spins configuration of Fig. ^] Each coupling is of dipole-dipole type, with the strength of the AB and 
BC couplings equal to 1 while for the AC coupling it is 1/8. The lack of any chemical shift makes the three spins 
indistinguishable ("full" controllability is therefore lost). The task is to have the three spins behave like a linear chain, 
canceling the effect of the AC coupling: 

^^fd = ^ (Alio + A220 — 2A330) 
— (Aoii + A022 — 2A033) 

Hs^-^ (Aioi + A202 - 2A303) 

by means of the control Hamiltonian 

He = u(Aooi + Aoio + Aioo). 

Since the purposes is now only to attain a useful control signal u, we can choose the initial condition ^(0) = ^d- 
Normally in the laboratory one uses g{0) — Pq(0) ® QpiO) ^ Q-yi^) where each reduced density is ahgned along the Ai 
axis. Using the feedback scheme of (|34() the so-called FID (free induction decay) signal + g''*^'^ + with and 
without control is shown in Fig. El The same tracking behavior holds for the entire tensor g. Hence the Hg coupling 




FIG. 12: Left: Uncontrolled FID response. Right: FID response with the feedback control. 



is effectively suppressed on the entire Hilbert space of the system. Due to the symmetries of the problem, in this 
case the local precession are totally ininfluent. No pulse sequence is known in the NMR literature able to achieve this 
decoupling, not even if we restrict to the one-dimensional submanifold corresponding to the FID signal and disregard 
what happens on the rest of the Hilbert space. This shows the potential of model-based feeback methods even for the 
usual task of open-loop control design. 



VII. CONCLUSION 



In the classical world, feedback can provide a degree of robustness to all tasks requiring state manipulation or 
protection which is by far unachievable by open loop control methods. If one adds that the complexity of the feedback 
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synthesis is much lower than that of open loop design [2g, it is easy to realize which one control engineers like it better. 
The main purpose of this paper is to show that to some extent known feedback methods could be used also in NMR 
systems. Unlike other contexts, where "quantum feedback" is studied 0,0,13123' complete noninvasivity of 
the bulk measurements makes it theoretically sound -at least in principle- to study the feedback problem for NMR 
systems in a purely deterministic and unitary framework. A possible feedback synthesis follows directly from standard 
control-theoretic Lyapunov techniques. 
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